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Capon’s  maximum  likelihood  method  (MLM)  is  being  used  increasingly  for  improved 
localization  and  high  resolution  of  sources  in  matched-field  processing  (MFP)  in  a  waveguide 
environment.  When  the  noise  component  is  dominated  by  modal  noise  (that  is,  noise  due  to 
excitation  of  the  modal  structure  of  the  waveguide  by  distant  sources  of  acoustic  energy,  e.g., 
wave  action,  distant  shipping),  the  MLM  can  become  unstable;  when  the  number  of  modes 
( M )  supported  by  the  waveguide  is  considerably  less  than  the  number  of  sensors  (TV),  a 
“reduced”  MLM  is  available  to  stabilize  the  processing  [Byrne  et  al.,  J.  Acoust.  Soc.  Am.  87, 
2493-2502  ( 1990)  ].  In  this  paper  the  sector-focused  stability  (SFS)  methods  of  Byrne  and 
Steele  (Proc.  IEEE.,  ICASSP  1987)  are  employed  to  treat  the  case  of  M nearly  equal  to  TV. 
Simulations,  using  a  shallow-water  range-independent  environment,  are  presented  to  illustrate 
the  increase  in  stability  of  SFS  over  MLM  in  the  presence  of  phase  errors  on  the  sensors  or 
more  general  mismatch. 

PACS  numbers:  43.30. Wi,  43.60.Gk 


INTRODUCTION 

Matched-field  processing  (MFP)  is  a  type  of  signal  pro¬ 
cessing  procedure  that  involves  the  comparison  of  signal  re¬ 
plica  vectors  obtained  from  models  of  the  acoustic  propaga¬ 
tion,  such  as  normal  mode  expansions  for 
range-independent  waveguides,  with  received  signal  data.  In 
cases  where  a  plane-wave  propagation  model  is  generally 
inadequate,  MFP  has  been  shown  to  be  a  useful  method  for 
detection  and  localization  of  sources.  The  need  to  go  beyond 
the  plane-wave  assumption  was  discussed  more  than  a  dec¬ 
ade  ago  by  Hinich1  and  Bucker,2  and  later  by  Mermoz.3 
Tindle  et  al.*  introduced  what  has  come  to  be  called  “mode 
space  processing.”  Hinich3  considered  modal  noise  as  a  su¬ 
perposition  of  sources  in  the  waveguide  and  showed  that  the 
statistical  maximum  likelihood  estimation  of  unknown  pa¬ 
rameters  could  be  achieved  with  the  number  of  sensors  (TV) 
equal  to  the  number  of  modes  (At)  supported  by  the  wave¬ 
guide.  Heitmcyer  et  al.6  considered  MFP  for  the  Pekeris 
model. 

To  achieve  high  resolution,  while  controlling  computa¬ 
tional  load,  several  authors  have  considered  nonlinear  meth¬ 
ods.  The  most  widely  studied  is  Capon’s  estimator,7  called 
the  “maximum  likelihood  method”  (MLM)  by  Lacoss,8  as 
well  as  the  “minimum  variance  distortionless  look”  by  oth¬ 
ers.9  For  a  detailed  discussion  of  MLM  applied  to  MFP  see 
Baggeroer  et  al.'°  Klemm"12  considered  a  maximum  en¬ 
tropy-type  estimator  and  studied  the  effects  of  mismatch  due 
to  incomplete  or  inaccurate  knowledge  of  environmental  pa¬ 


rameters,  such  as  those  describing  the  bottom.  Other  studies 
of  sensitivity  of  MFP  to  mismatch  have  considered:  pertur¬ 
bations  in  water  depth  and  geoacoustic  parameters;13  the 
use  of  waterborne  modes  only,  as  well  as  sound-speed  profile 
and  sediment  mismatch;14  and  reduction  of  array  aperture 
and  coverage  of  the  water  column.15  Several  authors  have 
investigated  the  estimation  procedures  themselves  and  have 
suggested  new  estimators  to  overcome  problems  with  mis¬ 
match.  Shang16-18  has  studied  the  mode  space  processing 
mentioned  in  Tindle  et  al.,*  and  Yang19  has  suggested  pro¬ 
jection  onto  certain  eigenvector  subspaces  to  increase  the 
number  of  resolvable  modes.  Shang20  has  also  discussed  the 
use  of  MFP  with  known  sources  to  determine  waveguide 
parameters. 

It  is  well  known  that  Capon’s  MLM  estimator  is  sensi¬ 
tive  to  mismatch  between  the  assumed  model  and  the  actual 
data.  For  the  plane-wave,  case  Byrne  and  Steele21  showed 
that  certain  eigenvectors  of  the  cross  spectral  matrix  (CSM) 
R  could  become  unstable  when  a  correlated  noise  compo¬ 
nent  is  present,  which  leads  to  deterioration  of  the  MLM  in 
the  presence  of  phase  errors  or  other  mismatch.  A  “sector- 
focused  stability”  (SFS)  method  has  been  developed22  to 
overcome  this  instability.  The  situation  in  the  normal  mode, 
rather  than  the  plane-wave,  case  is  closely  analogous  to  that 
studied  by  Byrne  and  Steele.21,22  The  presence  of  modal 
noise  (distant  noise  sources  exciting  the  modal  structure) 
creates  similar  eigenvector  instability.  Byrne  et  al.™  derived 
mode  space  MLM,  or  “reduced  ML”  (RML),  as  a  special 
case  of  SFS  and  showed  that  it  could  improve  stability  when 
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the  number  of  modes  is  significantly  less  than  the  number  of 
phones.  In  this  paper  we  develop  sector  focusing  further,  to 
achieve  increased  stability  and  reduced  computational  load 
in  the  normal  mode  case.  This  enhanced  sector  to  agree  with 
focusing  can  improve  RML  when  the  number  of  modes  is 
about  equal  to  the  number  of  phones. 

In  Sec.  I  we  present  some  background  on  the  issue  of 
sensitivity  to  mismatch  of  ML  and  on  the  recently  popular 
technique  of  dimensionality  reduction,  of  which  SFS  is  a 
special  case.  In  Sec.  II  we  present  the  SFS  method  in  detail, 
discuss  the  normal  mode  model,  and  consider  the  role  of 
normalization  in  matched-field  processing.  In  Sec.  Ill  we 
discuss  the  simulations  performed  and  describe  the  results. 
Finally,  we  have  a  brief  statement  of  the  conclusions  we 
draw  from  this  work. 

I.  BACKGROUND 

Shortly  after  the  appearance  of  Capon's  paper,7  Seleg- 
son24  showed  that  the  resolving  properties  of  the  MLM  esti¬ 
mator  could  deteriorate  if  there  are  phase  or  amplitude  er¬ 
rors  in  the  data.  He  also  noted  (in  a  private  communication 
recalled  by  McDonough25 )  that  the  wider  the  spread  be¬ 
tween  the  largest  and  smallest  eigenvalues  of  the  cross  spec¬ 
tral  matrix  R,  the  more  severe  the  deterioration.  McDon¬ 
ough25  obtained  quantitative  measures  of  the  sensitivity  of 
MLM  to  mismatch.  Cox26-27  viewed  sensitivity  to  mismatch 
as  related  to  ability  to  resolve  sources  or  to  reject  interfer¬ 
ence,  and  suggested  several  methods  to  improve  stability, 
including  multiple  main  lobe  constraints.  Vural2*  pointed 
out  that  mismatch  is  most  serious  at  high  signal-to-noise 
ratio  (SNR)  and  considered  the  use  of  multiple  main  lobe 
constraints  and  derivative  constraints  to  combat  the  prob¬ 
lem.  The  use  of  multiple  linear  constraints  was  also  treated 
by  Steele.2’  Gray50  noted  that  stability  can  be  obtained  by 
processing  the  output  of  several  presteered  beams  (beam 
space  processing),  provided  the  beams  are  sufficiently  or¬ 
thogonal.  More  recently,  Cox  et  al.3'  have  achieved  robust 
beamforming  by  imposing  a  constraint  on  white  noise  re¬ 
sponse. 

Byrne  and  Steele21  related  the  sensitivity  of  the  MLM 
estimator  to  that  of  particular  eigenvectors  of  the  CSM  and 
showed  that,  by  modifying  MLM  to  reduce  the  influence  of 
the  unstable  eigenvectors,  stable  estimation  was  possible. 
Their  analysis  showed  that  instability  increases  with  the 
presence  of  noise  components  that  appeared  signal-like  to 
the  array  (such  as  spherical  isotropic  noise  to  an  oversam¬ 
pled  line  array).  Of  course,  such  components  lead  to  the 
eigenvalue  spread  that  concerned  Seligson  and  McDon¬ 
ough.  To  combat  sensitivity,  Byrne  and  Steele22  suggested 
the  SFS  procedure,  which  also  has  the  advantage  of  dimen¬ 
sionality  reduction,  reducing  the  size  of  the  matrix  to  be  in¬ 
verted  for  estimation  by  using  an  N by  K(K<N)  projector 
matrix  V,  with  Vy  V  —  I  (the  f  denotes  conjugate  trans¬ 
pose)  ,  to  transform  R  to  the  K  by  K  matrix  Vf  R  V.  Dimen¬ 
sionality  reduction  has  become  popular  recently  and  has 
been  used  for  a  variety  of  purposes,  including  simple  reduc¬ 
tion  of  computational  load. 

Several  authors  have  considered  dimensionality  reduc¬ 
tion  for  the  plane-wave  case.  Bienvenu  and  Kopp22  suggest¬ 


ed  beam  space  processing  to  avoid  the  effects  of  nonwhite 
noise  on  eigenvector-based  methods.  Lee  and  Wengrovitz55 
showed  that  the  SFS  procedure  can  be  effective  in  overcom¬ 
ing  mismatch  due  to  finite  averaging  in  estimating  the  CSM. 
Van  Veen  and  Williams54  considered  general  dimensionality 
reduction  methods  for  reducing  computation,  and  Xu  and 
Buckley55  analyzed  the  statistical  improvement  due  to  di¬ 
mensionality  reduction.  Cox  et  al.36  have  studied  conven¬ 
tional  beamforming  on  selected  subarrays  for  MFP  with 
long  arrays,  to  reduce  computational  load. 

In  matched-field  processing  one  is  usually  dealing  with 
a  waveguide  that  imposes  structure  on  the  propagating  field. 
Distant  sources  of  acoustic  energy  (shipping,  surface 
waves)  can  excite  the  modal  structure  of  the  waveguide  and 
present  to  the  array  a  noise  vector  whose  structure  reflects 
that  waveguide.  We  refer  to  this  as  modal  noise.  When  the 
waveguide  supports  fewer  modes  than  phones  (M<N), 
such  modal  noise  presents  to  the  array  the  same  sort  of  struc¬ 
tured,  low-dimensional  vectors  presented  by  signal  sources. 
The  problem  is  closely  analogous  to  the  spherical  isotropic 
noise  and  oversampled  array  case  discussed  above.22 

Dimensionality  reduction  can  be  implemented  in  MFP 
through  so-called  mode  space  processing  (MSP).  In  a  relat¬ 
ed  paper,  Byrne  et  al.23  used  the  SFS  approach  to  derive 
mode  space  MLM  processing  and  to  examine  the  improve¬ 
ment  in  stability  this  procedure  affords.  When  the  number  of 
modes  is  not  significantly  less  than  the  number  of  sensors  the 
advantage  of  MSP  over  MLM  is  not  dramatic,  and  further 
reduction  of  dimension  is  warranted.  We  consider  a  method 
for  achieving  further  dimensional  reduction  in  this  paper: 
using  regions  of  range-depth  space  to  define  the  sectors;  and 
using  the  eigenstructure  of  related  positive  definite  matrices 
to  determine  the  essential  dimensions  of  the  sectors,  and  to 
provide  the  transformations  required. 

The  basic  idea  of  dimensionality  reduction  is  the  use  of  a 
nonsquare  N  by  K  matrix  V  (K<N)  with  orthonormal  co¬ 
lumns  ( Vy  V  =  /)  to  project  the  data  vectors  onto  a  subspace 
of  dimension  K.  The  resulting  CSM  of  the  projected  data  is 
T=  VfRV,  which  is  K  by  K.  All  subsequent  processing  in¬ 
volves  only  T.  If  the  columns  of  V  are  the  steering  vectors 
associated  with  pre-selected  beams,  then  Fis  the  CSM  used 
in  beam  space  processing.50  If  the  beams  are  not  orthogonal 
initially,  an  orthogonalization  procedure  (using,  say,  “QR 
factorization”  or  Gram-Schmidt57 )  is  needed  to  produce  V. 
If  the  columns  of  V  are  selected  columns  of  the  TVby  N  identi¬ 
ty  matrix,  then  F  is  the  CSM  obtained  by  deleting  certain 
sensors  of  the  array  (subarray  processing).  One  objective  of 
dimensionality  reduction  is  to  reduce  the  size  of  the  matrix 
to  be  inverted  without  destroying  much  signal  information. 
If  the  signal  vector  of  interest  lies  in  the  span  of  the  columns 
of  V,  then  there  is  no  loss  of  signal  energy.  Another  objective 
is  to  prewhiten  the  noise-only  CSM.  If  R  =  S  +  Q,  with  S 
and  Q  the  signals-only  and  noise-only  CSM,  respectively, 
then  to  whiten  Q  with  a  square  matrix  (no  dimensionality 
reduction )  requires  exact  knowledge  of  Q.  However,  VfQV 
can  be  made  close  to  a  K  by  K  identity  matrix,  when  K<N 
and  the  columns  of  Fare  obtained  from  nearby  beam-steer¬ 
ing  vectors,  that  is,  we  can  achieve  a  local  whitening  of  the 
noise  without  knowing  Q33 
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When  the  objective  is  to  stabilize  MLM  estimation,  the 
projection  matrix  V  is  selected  to  eliminate  or  reduce  the 
contribution  coming  from  unstable  eigenvectors  of  R.  For 
the  case  of  normal  mode  propagation  considered  in  this  pa¬ 
per,  the  columns  of  V  can  be  chosen  to  be  the  M  mode¬ 
sampling  vectors  (leading  to  the  RML  estimator23 )  or  can 
correspond  to  a  subsector  of  range  and  depth  space. 

II.  THE  NORMAL  MODE  MODEL,  SECTOR-FOCUSED 
STABILITY  AND  NORMALIZATION 

A.  The  normal  mode  model 

Using  the  notation  of  Shang16  we  write  the  field  excited 
by  a  point  source  in  the  far  field  of  a  waveguide  in  a  normal 
mode  expansion  as 

p(rjyZ0 )  =  iri  jr  Um(z)sm(rj0),  (1) 

m  —  1 

where  z0  is  the  source  depth,  z  is  the  receiver  depth,  and  r  is 
the  source  range  from  an  arbitrary  origin  (usually  where  the 
vertical  array  is  located).  The  U,„  ( z )  are  the  eigenfunctions 
of  the  depth-only  Sturm-Liouville  boundary  value  problem: 
Um(z)  is  called  the  rath  mode.  The  coefficients  s„(r^0) 
represent  the  excitation  of  the  mth  mode  by  a  source  at  range 
r  and  depth  z0 : 

s,„(rj0)  =exp(3ra'/4)exp(  -  P„r)  exp  (ikmr) 

X  Um  (zn  )^2rr/kmr ,  (2) 

where  km  and  /?„  are  horizontal  and  vertical  wave  numbers, 
respectively.  Higher  modes  and  continuous  modes  are  con¬ 
sidered  as  noise,  since  we  are  interested  in  far-field  sources. 
Sampling  the  pressure  field  using  sensors  at  depths  z„ 
(n  =  1 . .v;  we  obtain  the  (narrow-band)  data  vector 

P (r*z0  )  =  [p(r\z,  ),...,p(r,zNj„  )  ]  r.  (3) 

Letting  U  be  the  N  by  M  matrix  whose  rath  column  is 
[  ( U,„  (z,  ( z„ )  ]  r,  and  s (r,z0 )  be  the  M  by  1  vector 

whose  rath  entry  is  sm  (r,z0 ),  we  can  then  write 

p(r,z0)  =  Us(rfz0).  (4) 

We  assume  that  the  single  snapshot  data  vector  x  has  the 
form 

J 

x  =  ^  ^iP(rj<2oJ)  +  noise 

i  ■■  i 

=  2  Aiu,j  +  ur  +  v>  (5) 

J  -  l 

where  rJ  and  z0J  are  the  range  and  depth,  respectively,  of  the 
Jtb  source;  s,  =  s(rJfz0J);  the  entries  of  y  represent  the  exci¬ 
tation  of  the  various  modal  amplitudes  by  aggregate  noise 
sources;  rj  is  random  white  noise;  and  th eAy  are  independent 
random  variables  representing  time-varying  channel  fluctu¬ 
ations  (a  Rayleigh  fading  channel,  for  example).  To  lessen 
the  effects  of  noise  it  is  common  to  take  R  =  average  ( xxf ), 
where  the  averaging  is  performed  over  the  available  snap¬ 
shots.  Then  our  model  for  R  becomes 

R  =  USAS'U  t  +  l/Gt/+  +  H,  (6) 

where  the  jth  column  of  S  is  t(r,,zaj);  A  is  the  variance- 
covariance  matrix  of  the  A,  ;  UGU +  is  the  CSM  of  the  modal 
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noise  (so  G  —  average  (yy* )  and  has  the  same  form  as 
SAS* ,  but  for  a  large  J representing  the  contributions  of 
many  distant  noise  sources);  and  H  =  average  )  is  the 
CSM  of  nonmodal  noise. 

When  the  number  of  modes  M  is  less  than  the  number  of 
sensors  N,  both  the  signal  vectors  Usj  and  the  modal  noise 
vectors  Vy  represent  energy  projected  onto  a  smaller  dimen¬ 
sional  subspace,  the  span  of  the  columns  of  U,  so  that  signals 
and  modal  noise  are  not  easily  separated.  When  the  modal 
noise  component  UGU +  dominates  R,  the  lowest  eigenvec¬ 
tors  of  R  tend  to  be  nearly  orthogonal  to  all  the  columns  of  U 
and  hence  to  every  p  ( r\z0 )  =  Cs  ( r,z0 ).  This  is  closely  analo¬ 
gous  to  the  plane-wave  situation,  with  spherical  isotropic 
noise  and  an  oversampled  linear  array,  and  will  cause  MLM 
to  be  sensitive  to  phase  errors. 

If  M  is  much  less  than  N,  stability  can  be  achieved  by 
letting  the  V  in  the  SFS  method  be  the  (orthogonalized) 
columns  of  U.2X  If  M  is  not  much  less  than  N,  greater  stabil¬ 
ity  can  be  obtained  by  choosing  as  our  sector  a  subset  2  of 
0  =  (r,z)  space  whose  essential  dimension  is  K  <M,  and 
then  proceeding  with  SFS  on  the  reduced  K  by  K  matrix 
yf  R  V,  as  described  below. 

In  Sec.  Ill  we  discuss  simulated  cases  of  normal  mode 
propagation  and  the  use  of  SFS  to  stabilize  the  MLM  estima¬ 
tor. 

B.  Sector-focused  stability 

We  assume,  initially,  an  arbitrarily  configured  array  of 
N  sensors.  The  narrow-band  single-snapshot  data  vector  is 
the  N by  1  vector  x  and  the  (sampled)  cross  spectral  matrix 
is  R  —  average  (xx+  ),  where  the  average  is  taken  over  the 
available  number  of  snapshots.  We  assume  that  R  has  (ap¬ 
proximately)  the  form 

R=  £  +  Q,  (7) 

i=  i 

where  the  p(0( )  are  the  signal  vectors  assumed  to  lie  in  a 
parametrized  set  {p(0),  6  in  fl}  and  Q  is  the  noise-only 
CSM.  For  now  the  0  have  no  further  significance  and  maybe 
scalar  or  vector  quantities.  The  objective  is  to  determine  J ,  as 
well  as  the  6 )  and  (perhaps)  the  a2,  usually  without  knowl¬ 
edge  of  Q.  Capon’s  MLM  estimator  is  a  popular  choice  be¬ 
cause  it  does  not  require  Q  and,  for  sufficiently  high  signal- 
to-noise  ratio  (SNR),  achieves  better  resolution 
(determination  of  J)  than  conventional  methods  (e.g.,  Bart¬ 
lett).7’*  Other  procedures,  such  as  “optimal  array  process¬ 
ing"  or  eigenvector-based  methods,  require  a  prewhitening 
of  R  that  necessitates  knowledge  of  Q. 

Capon’s  MLM  proceeds  as  follows:  for  each  fixed  value 
0  =  a)  we  find  the  linear  filter  f  =  f(<u)  that  minimizes  aver¬ 
age  output  power  f  R  f,  subject  only  to  the  constraint  of  no 
distortion  of  p(o>)  =  1.  The  well-known  solution  is 

f(< o)=R  lp(<u)/[p(ra)tR  ’plat)].  (8) 

Applying  this  optimal  filter  to  the  actual  data  gives  the 
average  output  power  f*  (w)  R  f{co)  equal  to 

l/[p  (w)7/?  'p(<u)]  =  MLM(w),  (9) 

which  is  Capon's  estimate  of  the  power  associated  with 
9  =  a.  Because  of  the  appearance  of  R1  in  (9),  those  ei- 
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gen  vectors  of  R  corresponding  to  the  lowest  eigenvalues  ( let 
us  call  them  the  “lowest”  eigenvectors)  have  the  greatest 
influence  on  MLM(&>).  As  discussed  by  Byrne  and  Steele,21 
the  presence  of  a  correlated  noise  component  Q,  correspond¬ 
ing  to  the  superposition  of  many  signal-like  noise  sources, 
can  cause  these  lowest  eigenvectors  to  become  unstable  in 
the  presence  of  mismatch  (such  as  phase  errors).  To  stabi¬ 
lize  the  estimation  of  the  0j  one  must  reduce  or  eliminate  the 
contribution  from  these  unstable  eigenvectors.  The  SFS 
method  is  designed  to  achieve  this  stability. 

The  MLM  estimator  is  the  averaged  output  of  the  linear 
filter  in  (8),  which  is  dominated  by  the  lowest  eigenvectors 
of  R.  When  the  modal  noise  ( UGW )  component  of  Q  is 
significant,  the  lowest  eigenvectors  of  R  tend  to  be  nearly 
orthogonal  to  the  columns  of  U,  hence  to  all p(0),  and  MLM 
becomes  unstable.  The  sector-focused  stability  (SFS)  meth¬ 
od  employs  a  modification  of  ( 8 )  to  eliminate  the  instability, 
at  least  for  estimation  within  a  pre-selected  subset  2  of  (l. 

For  any  Nby  K  matrix  V,  every  linear  filter  g  =  g(u>) 
can  be  written  uniquely  as  g  =  Va  +  w,  where  Vtw  =  0  and 
a  is  some  K  by  1  vector.  The  more  w  dominates  g  the  more  g 
will  be  orthogonal  to  the  columns  of  V.  Those  g  whose  w 
component  is  missing  are  those  least  orthogonal  to  the  co¬ 
lumns  of  V;  thus  the  additional  constraint  that  g(o>)  =  Va 
for  some  a  simply  makes  it  harder  for  g(<u)  to  be  orthogonal 
to  the  columns  of  V.  The  answer  is  to  build  V  from  vectors 
p(0),  0  in  2,  thus  making  g(co)  more  stable  for  estimation 
within  2. 

To  obtain  the  projection  matrix  V  we  begin  by  simulat¬ 
ing  a  noise  component  equal  to  the  superposition  of  many 
independent  sources  within  2;  that  is,  we  form  (an  approxi¬ 
mation  of)  Q(2)  —  average  {p(0)p(#)+},  for  0  in  2.  Tak¬ 
ing  the  eigenvector/eigenvalue  decomposition  Q(2) 
—  WL  fV\  W*  W  —  I,  L  =  diag  {A,  ••  -A„  },  we  define  the 
essential  dimension  of  2  to  be  the  number  k  of  (essentially) 
nonzero  terms  in  L.  The  matrix  V  has  for  its  columns  these  k 
columns  of  W,  so  V  is  N  by  K  and  V1  V  —  I. 

We  derive  the  SFS  method  by  modifying  Capon’s  meth¬ 
od  as  follows:  for  each  fixed  0  —  to  within  2  we  find  the  linear 
filter  g  =  g(ru)  that  minimizes  g*  R%  subject  to  the  con¬ 
straints  g+p (to)  =  1  and  g  =  Va,  where  a  is  a  AT  by  1  vector 
and  therefore  g  is  in  the  span  of  the  columns  of  V.  It  is  unim¬ 
portant  what  a  is;  the  point  is  that  g (to)  cannot  now  have  a 
sizable  component  orthogonal  to  the  columns  of  V.  Proceed¬ 
ing  as  before,  we  find  that  the  optimal  filter  is 

g(ru)  =  V(VfRV)  ~lVfp(m) 

X  lp(<u)fF(  VfR  V)  ~  1  Vfp(to)  ]  ~ 1  (10) 

and  the  average  output  power  is  the  SFS  estimator: 

t?(o>)Rg(6>)  =  l/[p {co)'V(V'RV)  ~  'V'plo)] 

=  SFS(fu).  (11) 

C.  Rota  of  normalization  In  matched-field  processing 

As  they  appear  in  Eqs.  (9)  and  (11),  MLM(<»)  and 
SFS(<y)  have  singularities  wherever  p(&)  approaches  zero. 
This  is  principally  a  problem  near  the  surface  of  a  waveguide 
since,  from  die  acoustic  boundary  conditions,  p(  to)  must 
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vanish  at  the  surface.  This  is  easily  seen  in  Eq.  ( 1 ).  There¬ 
fore,  if  one  wishes  peaks  in  the  ambiguity  function  to  reflect 
the  presence  of  acoustic  sources  and  not  the  relative  magni¬ 
tude  of  the  individual  p(o),  then  a  normalization  of  the  pro¬ 
cessor  is  required. 

Before  deciding  what  this  normalization  should  be,  we 
must  consider  the  processor  response  to  the  noise-only  case, 
to  see  what  background  the  processor  will  provide.  First, 
consider  the  case  of  white  noise,  which  is  totally  uncorrelat¬ 
ed  between  the  hydrophones  so  that  the  time-averaged  cross- 
spectral  matrix  R  =  Q  is  completely  diagonal.  Further,  as¬ 
sume  that  the  average  acoustic  power  delivered  to  each 
hydrophone  is  constant  so  that  R  is  some  multiple  of  the 
identity  matrix.  In  this  case  it  is  simple  to  evaluate  MLM  ( co ) 
from  (9): 

MLM(o)  =  W(«)/p+(<u)p(«).  ( 12) 

where  Af(o>)  is  the  (as  yet  undetermined)  normalization 
function.  The  sector-focused  estimator  SFS(&>)  (11)  be¬ 
comes 

SFS(w)  =N(o>)/pHa>)V(V*V)  -  'K+p(«>)>  (13) 

=  JV(*>)/pt(v)VVfP(v),  (14) 

since  K+  V  =  /. 

Before  interpreting  what  we  have  just  found,  let  us  con¬ 
sider  the  other  extreme  type  of  noise  field  one  might  encoun¬ 
ter  in  an  acoustic  waveguide.  Let  us  consider  noise  propagat¬ 
ing  to  the  array  from  a  large  number  of  distant  point  sources 
randomly  distributed  in  depth  and  range  (in  fact,  distribu¬ 
tion  in  range  only  is  sufficient)  so  that  no  near-field  effects 
are  included.  Then  it  has  been  shown  that  the  mode  space 
cross-spectral  matrix  will  be  diagonal  in  form.2-1'"*'40  If  we 
also  assume  that  the  average  acoustic  power  carried  by  each 
normal  mode  is  equivalent,  then  the  mode  space  CSM  [i.e., 
G  in  Eq.  ( 6)  ]  is  simply  a  multiple  of  the  identity  matrix  and 
the  phone  space  matrix  is  R  =  £?  =  UU* .  Because/?  =  UU* 
is  not  invertible,  we  cannot  resort  to  (8)  and  (9)  for  the 
MLM.  We  must  directly  solve  the  problem:  minimize 
f*  (.(o)R  f(<u)  =  f*  (eo)  l/l/+f(zt>), subject  tof*  (<y)p(ru)  =  1. 
The  MLM  esimate  will  then  be  f*(u))/f  f(<w)  for  the  optimal 
f(«).  Writingg  =  U*  f(<u),  the  problem  becomes:  minimize 
g^subjectto  1  =  f*(<u)p(<y)  =  f  (a)  Us(a)  =  g+s(«).The 
solution  is 

g=  [s^wlsfu))]  's(<a).  (15) 

From 

g'g  =  f(a>)RHa>)  =f(o>)UU'f(a>),  (16) 

we  see  that  the  MLM  is  gf  g  or, 

MLM(<u)  =  Ar(<u)/st(<w)s(a>).  (17) 

Similarly,  SFS(o>)  takes  the  form, 

SFS(«)  =  N(a)/a\<a)U'V{VWU'V)  ~'VfVsM.  (18) 

If  the  columns  of  Van  simply  orthogonalized  columns 
of  U  then  V  =  UG  for  some  invertible  G,  and  it  follows  that 
SFS(<u)  reduces  to 

SFS(<u)  =  N(o)W(<o)s(o>).  ( 19) 

Now  that  we  know  what  MLM(<o)  and  SFS(<u)  look  like  in 
these  noise  fields,  we  turn  to  the  question  of  the  appearance 
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of  their  corresponding  processor  responses. 

For  the  ambiguity  surface  to  be  useful,  the  background 
against  which  signals  are  sought  must  be  relatively  flat.  If  the 
noise  component  is  white,  with  the  MLM  surface  for  white 
noise  only  being  given  by  (12),  the  best  normalization  is 
N(co)  —  p(co)f  p(o}).  However,  if  the  noise  is  mainly  modal, 
normalizing  MLM  by  N(<o)  =  st(<u)s(o>)  will  produce  a 
flatter  background.  Because  we  are  concerned  in  this  paper 
with  the  effects  of  a  sizable  modal  noise  component  we  have 
chosen  the  latter  normalization. 

III.  SIMULATIONS 
A.  Environmental  model 

In  order  to  demonstrate  the  effectiveness  of  the  sector- 
focused  stability  technique,  we  have  performed  numerical 
simulations  of  a  shallow-water,  range-independent  environ¬ 
ment,  as  shown  schematically  in  Fig.  1.  This  environment, 
which  is  modeled  after  a  location  on  the  continental  shelf  off 
the  southern  California  coast,  consists  of  a  water  channel  of 
500-m  depth  overlying  a  sediment  layer  and  an  isospeed, 
semi-infinite  subbottom.  The  sound-speed  profile  for  the  wa¬ 
ter  and  sediment  layers  and  top  section  of  the  subbottom  are 
shown  in  Fig.  2.  In  addition  to  the  sound-speed  difference, 
there  is  also  a  density  contrast  between  the  water  layer  and 
the  sediment.  The  water  density  is  1 .0  g/cc  and  the  sediment 
density  is  1 .5  g/cc.  There  is  an  attenuation  factor  in  the  sedi¬ 
ment  layer  of  0. 1  dB/A. 

The  acoustic  field  calculations  for  these  simulations 
were  performed  using  the  SACLANTCEN  normal-mode 
acoustic  propagation  (SNAP)  model.  For  the  set  of  environ¬ 
mental  parameters  shown  in  Figs.  1  and  2,  the  SNAP  model 
indicates  14  propagating  modes  at  a  source  frequency  of  30 


PIO.  I .  Schematic  diagram  showing  the  acoustic  environmental  model  used 
in  the  simulations.  The  geometrical  arrangement  of  the  array,  the  place¬ 
ment  of  surface  noise  sources,  and  the  delineation  of  the  search  space  (later 
divided  into  search  sectors)  are  also  indicated  in  this  diagram. 


Hz.  The  vertical  receiving  array  consists  of  3 1  equally  spaced 
hydrophones  spanning  the  central  450  m  of  the  water  col¬ 
umn.  Here  we  will  be  concerned  with  propagation  over 
ranges  between  10  and  1 8  km  from  the  source. 


B.  Construction  of  the  CSM 

The  matched-field  processor  accepts  as  data  input  the 
cross  spectral  matrix  R,  which  is  calculated  as  the  sum  of  the 
signal-only  matrix  S  and  the  noise-only  matrix  Q, 

R=S+Q.  (20) 

For  a  source  located  at  range  r  and  depth  z0,  SNAP  yields 
p(r,z0 )  and  s(r,z0 )  of  Eq.  (9),  from  which  S  may  be  calcu¬ 
lated, 

5  =  p(/vz0)p+(/vz0).  (21) 

The  magnitude  of  the  elements  of  p(r0,z0 )  are  adjusted  to 
scale  the  trace  of  this  matrix  to  give  the  signal  level  required 
for  the  simulation. 

The  noise  matrix  Q  has  two  components.  The  first  is  the 
“modal  noise”  matrix,  representing  noise  that  is  correlated 
between  the  hydrophones  when  detected  at  the  array.  This  is 
modeled  by  summing  contributions  from  1500  individual 


SOUND  SPEED  (m/s) 


FIG.  2.  Sound-speed  profile  for  the  continental  shelf  environment  used  in 
the  simulations. 
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sources,  which  are  placed  at  a  depth  of  7  m  and  distributed 
evenly  in  ranges  from  0.3  to  500  km  from  the  array.  The 
SNAP  model  implicitly  assumes  cylindrical  symmetry  when 
calculating  the  acoustic  field.  Therefore  there  is  a  require¬ 
ment  that  each  of  the  1500  sources  in  the  sum  represent  an 
annulus  of  sources  surrounding  the  array  whose  radius  is  the 
range  of  the  source  from  the  array.  To  accomplish  this,  the 
input  level  of  the  sources  must  be  scaled  so  that  the  average 
intensity  registered  at  the  hydrophones  from  all  of  them  will 
be  the  same.  The  second  component  of  Q  is  noise  that  is 
completely  uncorrelated  between  the  hydrophones  (“white 
noise”)  and  may  be  modeled  as  a  multiple  p  of  the  identity 
matrix.  The  total  noise  matrix  then  becomes 

1500 

e=£p,p  l+pl-  (22) 

1=  1 

In  these  simulations  we  will  monitor  the  relative  degradation 
of  the  MLM  and  SFS  estimators  when  phase  errors  are  intro¬ 
duced  onto  the  signals  received  at  the  hydrophones  array. 
The  introduction  of  these  phase  errors  takes  the  form  of  a 
similarity  transformation  of  R, 

R=DRD\  (23) 

where  D  =  diag  {e'"'},  and  the  8J  are  derived  from  a  set  of 
random  numbers  which  are  scaled  to  give  the  desired  maxi¬ 
mum  phase  error  incident  upon  the  array.  Figure  3  shows  a 
realization  of  the  random  set  0},  scaled  to  give  a  maximum 
value  for  the  error  of  ±  1  deg.  From  Eq.  (23)  we  see  that  the 
similarity  transformation  amounts  to  a  rotation  of  the  data 
matrix  about  a  direction  determined  by  the  relative  magni¬ 
tude  of  the  phase  errors  and  through  an  angle  determined  by 
their  overall  magnitude.  Although  the  perturbations  consid¬ 
ered  here  were  phase  errors,  other  forms  of  mismatch  would 
also  lead  to  similar  degradation  of  MLM.  The  SFS  approach 
offers  greater  stability  to  other  types  of  mismatch,  such  as 
erroneous  environmental  parameters,  since  the  exact  nature 
of  the  perturbations  is  not  important  in  the  development  of 
SFS. 

In  this  paper  we  use  four  parameters  to  characterize  the 
simulated  data  cross-spectral  matrix.  First,  the  source  level 
(SL)  is  the  average  of  the  diagonal  elements  of  the  signal- 
only  matrix  S,  expressed  in  dB  relative  to  unity.  Second,  the 
total  noise  level  (TNL)  is  defined  as  the  average  diagonal 
element  of  the  total  noise  matrix  Q,  also  expressed  in  dB 
relative  to  unity.  Third,  the  modal  to  white  noise  ratio 
( M  WR )  is  the  ratio  of  modal  to  white  noise  power,  again  in 
dB.  Finally,  the  random  phase  error  (RPE)  identifies  the 
magnitude  of  the  phase  errors,  and  is  expressed  as  a  multiple 
of  the  random  number  elements  displayed  in  Fig.  3.  For  ex¬ 
ample,  an  RPE  =  10  deg  means  that  we  multiply  the  phase 
errors  in  Fig.  3  by  10,  which  would  give  values  of  dj  varying 
between  —  lOdegand  +  lOdeg,  and  then  form  the  rotation 
matrix  D. 

C.  Comparison  of  MLM  and  SFS 

Now  that  we  have  the  simulated  data  matrix  R,  the 
problem  is  to  attempt  to  extract  the  source  information  from 
it  using  matched-field  techniques.  Here  we  will  compare  the 
performance  of  Capon’s  MLM  to  SFS  when  phase  errors  are 


FIG.  3.  Realization  of  the  set  of  random  phases  (maximum  value  =  ±  1 
deg)  used  to  introduce  phase  errors  onto  the  hydrophone  array  by  means  of 
the  similarity  transform  (23). 


present  in  the  data.  Since  we  are  concerned  only  with  stabi¬ 
lizing  Capon’s  MLM  we  do  not  include  comparisons  with 
conventional  ( Bartlett )  estimation.  The  advantages  and  dis¬ 
advantages  of  MLM  and  Bartlett  have  been  discussed  at 
length  in  the  literature.8,10 

An  important  input  to  all  matched-field  processors  is 
the  set  of  replica  field  vectors  (or  steering  vectors)  p(r,  ,z, ), 
where  /  denotes  the  /th  search  point  (or  trial  source  loca¬ 
tion).  Again,  we  have  used  SNAP  to  produce  these  vectors 
for  a  grid  of  points  within  a  search  space  selected  to  span  0  to 
200  m  in  depth  and  lOto  18min  range  (see  Fig.  1).  Now  for 
each  range/depth  point  in  the  search  grid  we  can  calculate 
MLM  using  Eq.  (9): 


MLM(r,,z, ) 


sV,,z,)s(r,,z,) 
p+(r1,z,)/{ ’p(r,,2,)  ’ 


(24) 


where  we  have  implemented  the  sfs  normalization  of  Eq. 

(17). 

When  calculating  the  SFS  estimator,  one  is  free  to 
choose  the  size  and  shapes  of  the  individual  sectors.  We  leave 
the  question  of  optimal  sector  geometry  for  future  investiga¬ 
tion.  Our  aim  here  is  to  show  that  even  a  very  simple  choice 
of  sectors  yields  a  very  stable  processing  result.  In  our  simu¬ 
lations  of  SFS  we  divided  the  total  search  space  into  16 
equally  sized  sectors.  Each  sector  covers  50  m  in  depth  and  2 
km  in  range.  Thus,  calculation  of  the  sector-focused  proces- 
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sor  requires  the  formation  of  16  matrices  V  describing  each 
of  the  16  sectors  within  the  search  space. 

For  each  sector  2,  the  calculation  of  V proceeds  in  three 
main  steps.  First  we  form  a  cross-spectral  matrix  as  follows: 


Q(2)=  X  (25) 

many  sources 

where  the  sum  is  taken  over  many  sources  distributed 
throughout  the  sector  of  interest.  The  exact  number  of 
sources  in  this  sum  is  not  important  so  long  as  it  is  large 
compared  to  the  number  of  hydrophones  in  the  receiving 
array.  This  assures  that  the  eigenvectors  of  N  will  span  the 
sector  in  question. 

Next,  we  numerically  determine  the  eigenvalues  and  ei¬ 
genvectors  of  g(2).  The  columns  of  V are  then  chosen  to  be 
a  set  of  K  eigenvectors  corresponding  to  the  largest  (not 
close  to  zero )  eigenvalues  of  Q{  2 ) .  The  exact  value  of  K  and 
the  form  of  V  depends,  in  general,  on  both  the  size  and  loca¬ 
tion  of  the  sector  in  question.  However,  K  is  limited  to  take  a 
value  between  one  and  the  number  of  propagating  modes  in 
the  waveguide.  This  is  easily  demonstrated  by  considering 
two  limiting  cases.  First,  consider  shrinking  the  size  of  the 
sector  to  a  single  point.  In  this  case  Q{2)  is  composed  of  a 
single  dyad  and  will  therefore  have  exactly  one  nonzero 
eigenvalue.  Second,  expand  the  sector  to  span  the  entire 
waveguide.  Here  the  number  of  large  eigenvectors  will  be 
equal  to  the  number  of  degrees  of  freedom.  For  an  acoustic 
waveguide  this  is  equal  to  the  number  of  modes,  and  so  the 
matrix  Q(  2 )  will  have  M  nonzero  eigenvalues.  In  our  simu¬ 
lations  we  choose  K  —  1 1,  so  that  V  is  a  31  X  1 1  matrix. 

Within  each  sector,  we  can  now  form  the  SFS  estimator 
as  follows: 


SFS(r,  ,z, ) 


_ st(r,,zl)s(r„z,) _ 

p \r„z,)V(v'R'V)  ‘ 


(26) 


D.  Quantification  of  the  estimator  functions 

In  order  to  quantify  the  performance  of  the  estimators, 
and  to  distinguish  between  visually  similar  ambiguity  sur¬ 
faces,  we  calculate  a  detection  factor  for  each  function  as 
follows: 

Detection  Factor:  PBR  =  (P  expressed  in  dB), 

(27) 

where  P  is  the  value  of  the  ambiguity  function  at  the  known 
source  location,  n  is  the  average  level  of  the  ambiguity  func¬ 
tion  ignoring  the  region  immediately  around  the  source 
peak,  and  a  is  the  standard  deviation  of  the  ambiguity  func¬ 
tion  ignoring  the  region  immediately  around  the  source 
peak. 

This  factor  has  been  discussed  previously  in  the  litera¬ 
ture.  1 1  We  do  not  infer  statistical  behavior  from  it,  but  use  it 
as  a  measure  of  “goodness.” 


E.  Results  and  analysis 

In  all  of  the  simulations  considered  here  we  maintain 
constant  values  for:  the  signal  level  (SL  =  50  dB);  the  total 
noise  level  (TNL  =  60  dB);  and  the  modal-to-white  noise 
ratio  (MWR  =  30  dB).  We  also  consider  a  single  source 
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FIG.  4.  MLM  ambiguity  surface  with  a  random  phase  error  of  0  deg. 


located  at  a  range  of  14.5  km  and  a  depth  of  30  m.  We  will  be 
concerned  with  the  performance  of  the  MLM  and  SFS  esti¬ 
mators  as  a  function  only  of  the  random  phase  error  intro¬ 
duced  at  the  hydrophone  array.  Essentially,  the  requirement 
here  is  to  localize  a  source  with  input  level  —  10  dB  within  a 
noise  field  dominated  by  spatially  correlated  ( modal )  noise. 

In  Fig.  4  we  see  the  MLM  ambiguity  surface  for  an 
RPE  =  0  deg.  On  this  surface  the  signal  peak  is  clearly  resol¬ 
vable  against  a  moderate  background  with  a  detection  factor 
PBR  =  9.06  dB,  and  is  located  at  the  correct  range  and 
depth.  This  is  to  be  expected,  since  there  is  no  difference 
between  the  replica  field  vectors  and  those  used  to  construct 
the  data  matrix.  In  Fig.  5  we  see  the  SFS  ambiguity  surface, 
also  for  -he  case  of  zero  phase  error.  This  surface  was  pro¬ 
duced  by  using  the  division  of  the  search  region  into  16  sec¬ 
tors,  each  50  m  in  depth  by  2  km  in  range,  described  above. 
This  figure  indicates  that,  in  the  RPE  =  0  deg  case,  the  SFS 
estimator  offers  no  advantage  over  MLM.  Although  the  sig¬ 
nal  peak  is  still  easily  resolvable  and  correctly  located,  the 
background  is  actually  slightly  worse  than  with  MLM.  The 
value  of  PBR  falls  to  8.28  dB. 

Let  us  now  see  the  effects  of  increasing  the  phase  error. 
In  Fig.  6  we  display  the  MLM  ambiguity  surface  for  an 
RPE  =  45  deg.  It  is  seen  that,  with  this  amount  of  phase 
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FIG.  6.  MLM  ambiguity  surface  with  a  random  phase  error  of  45  deg. 


error,  the  MLM  processor  has  been  degraded  and  distorted 
to  the  point  where  the  signal  peak  has  vanished  and  has  been 
replaced  by  two  large  spurious  sidelobes.  One  of  these  is  lo¬ 
cated  at  range  1 1.5  km  while  the  other,  although  found  at 
range  14.S  km,  is  located  at  the  surface  rather  than  at  the 
input  depth  of  30  m.  Evidently  the  MLM  processor  is  inef¬ 
fective  for  localization  purposes  once  this  level  of  phase  error 
is  reached.  In  contrast,  we  see  the  SFS  ambiguity  surface  for 
an  RPE  =  45  deg  in  Fig.  7.  Here  the  signal  peak  is  still  clear¬ 
ly  resolvable  above  the  background.  In  fact,  comparison  of 
Fig.  7  with  Fig.  5  shows  that  the  introduction  of  phase  error 
has  led  to  no  discernible  degradation  or  distortion  of  the  SFS 
surface.  Rather,  the  quality  of  the  surface,  from  a  signal  lo¬ 
calization  viewpoint,  seems  to  have  improved  slightly,  with 
the  value  of  PBR  increasing  to  8.55  dB  when  the  RPE  is  45 
deg. 

In  Fig.  8  curves  are  plotted  that  indicate  the  variation  of 
PBR  for  both  the  MLM  and  SFS  processors  as  a  function  of 
RPE.  As  the  phase  error  is  increased  we  see  the  MLM  curve 
fall  very  steeply  from  its  initial  value  of  9.06  dB  at  RPE  =  0 
deg,  to  a  value  of  2.55  dB  at  RPE  =  20  deg.  It  then  decreases 
more  slowly  as  the  RPE  is  further  increased,  reaching  1.81 
dB  at  RPE  =  100  deg.  However,  the  SFS  curve  is  dramati¬ 
cally  different.  As  the  phase  error  is  increased,  the  value  of 
PBR  also  increases  (see  again  Figs.  5  and  7),  reaching  its 
maximum  of  8.55  dB  at  RPE  =  45  deg.  After  this  point,  the 
values  of  PBR  begin  tq  slowly  fall;  but  the  processor  still 


FIG.  7.  SFS  ambiguity  surface  with  a  random  phase  error  of  43  deg. 
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FIG.  8.  Variation  of  PBR  for  both  the  MLM  and  SFS  processors  as  a  func¬ 
tion  of  random  phase  error. 


performs  effectively  out  to  relatively  large  values  of  RPE, 
giving  a  PBR  of  5.99  dB  at  RPE  =  90  deg  and  a  PBR  of  3.84 
dB  at  RPE  =  100  deg.  Apart  from  a  narrow  range  near 
RPE  =  0  deg  ( where  both  processors  work  almost  identical¬ 
ly)  the  SFS  consistently  outperforms  the  MLM  by  2-6  dB 
over  the  range  of  RPE  considered. 

IV.  CONCLUSIONS 

An  improved  matched-field  processor,  the  “sector-fo¬ 
cused  stability”  (SFS)  method,  has  been  developed  to  pro¬ 
vide  more  stable  source  estimation  performance  in  a  wave¬ 
guide  in  the  presence  of  high  levels  of  correlated  (modal) 
noise.  Under  these  conditions  the  “maximum  likelihood” 
(ML)  method  shows  serious  deterioration  whe  i  small  ran¬ 
dom  phase  errors  are  introduced  onto  the  hydrophones. 
However,  the  SFS  method  continues  to  return  stable  and 
accurate  source  estimates  with  levels  of  random  phase  error 
which  are  several  times  that  which  induces  failure  of  the  ML 
method. 
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